assignment 1

years <- 2017:2019
quarters <- 1:4
type= "Electric"

pge_elec <- NULL

for(year in years) {
  for(quarter in quarters) {
      filename <- 
        paste0(
          "PGE_",
          year,
          "_Q",
          quarter,
          "_",
          type,
          "UsageByZip.csv"
          )

      print(filename)
      temp <- read_csv(filename)
      pge_elec <- rbind(pge_elec,temp)
    }
}

## [1] "PGE_2017_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2017_Q2_ElectricUsageByZip.csv"
## [1] "PGE_2017_Q3_ElectricUsageByZip.csv"
## [1] "PGE_2017_Q4_ElectricUsageByZip.csv"
## [1] "PGE_2018_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2018_Q2_ElectricUsageByZip.csv"
## [1] "PGE_2018_Q3_ElectricUsageByZip.csv"
## [1] "PGE_2018_Q4_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q2_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q3_ElectricUsageByZip.csv"
## [1] "PGE_2019_Q4_ElectricUsageByZip.csv"

year = 2020
quarters <- 1:2
type= "Electric"
  for(quarter in quarters) {
      filename <- 
        paste0(
          "PGE_",
          year,
          "_Q",
          quarter,
          "_",
          type,
          "UsageByZip.csv"
          )

      print(filename)
      temp <- read_csv(filename)
      pge_elec <- rbind(pge_elec,temp)
  }

## [1] "PGE_2020_Q1_ElectricUsageByZip.csv"
## [1] "PGE_2020_Q2_ElectricUsageByZip.csv"

pge_elecfilter <- 
  filter(pge_elec, 
         CUSTOMERCLASS %in% 
           c("Elec- Residential","Elec- Commercial")) %>% 
  mutate(
    TOTALKBTU =
      TOTALKWH*3.4121
  ) %>% 
  select(
    !c(COMBINED, AVERAGEKWH, TOTALKWH)
  )

years <- 2017:2019
quarters <- 1:4
type= "Gas"

pge_gas <- NULL

for(year in years) {
  for(quarter in quarters) {
      filename <- 
        paste0(
          "PGE_",
          year,
          "_Q",
          quarter,
          "_",
          type,
          "UsageByZip.csv"
          )

      print(filename)
      temp <- read_csv(filename)
      pge_gas <- rbind(pge_gas,temp)
    }
}

## [1] "PGE_2017_Q1_GasUsageByZip.csv"
## [1] "PGE_2017_Q2_GasUsageByZip.csv"
## [1] "PGE_2017_Q3_GasUsageByZip.csv"
## [1] "PGE_2017_Q4_GasUsageByZip.csv"
## [1] "PGE_2018_Q1_GasUsageByZip.csv"
## [1] "PGE_2018_Q2_GasUsageByZip.csv"
## [1] "PGE_2018_Q3_GasUsageByZip.csv"
## [1] "PGE_2018_Q4_GasUsageByZip.csv"
## [1] "PGE_2019_Q1_GasUsageByZip.csv"
## [1] "PGE_2019_Q2_GasUsageByZip.csv"
## [1] "PGE_2019_Q3_GasUsageByZip.csv"
## [1] "PGE_2019_Q4_GasUsageByZip.csv"

year = 2020
quarters <- 1:2
type= "Gas"
  for(quarter in quarters) {
      filename <- 
        paste0(
          "PGE_",
          year,
          "_Q",
          quarter,
          "_",
          type,
          "UsageByZip.csv"
          )

      print(filename)
      temp <- read_csv(filename)
      pge_gas <- rbind(pge_gas,temp)
  }

## [1] "PGE_2020_Q1_GasUsageByZip.csv"
## [1] "PGE_2020_Q2_GasUsageByZip.csv"

pge_gasfilter <- 
  filter(pge_gas, 
         CUSTOMERCLASS %in% 
           c("Gas- Residential","Gas- Commercial")) %>% 
  mutate(
    TOTALKBTU =
      TOTALTHM * 100
  ) %>% 
  select(
    !c(COMBINED, AVERAGETHM, TOTALTHM)
  )

pge_consump <- rbind(pge_elecfilter,pge_gasfilter)

ca_counties <- counties("CA", cb = T, progress_bar = F)

bay_county_names <-
  c(
    "Alameda",
    "Contra Costa",
    "Marin",
    "Napa",
    "San Francisco",
    "San Mateo",
    "Santa Clara",
    "Solano",
    "Sonoma"
  )

bay_counties <-
  ca_counties %>%
  filter(NAME %in% bay_county_names)

usa_zips <- 
  zctas(cb = T, progress_bar = F)

bay_zips <-
  usa_zips %>% 
  st_centroid() %>% 
  .[bay_counties, ] %>% 
  st_set_geometry(NULL) %>% 
  left_join(usa_zips %>% select(GEOID10)) %>% 
  st_as_sf() 

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude data

bay_consump  <-
   pge_consump %>% 
   filter(ZIPCODE %in% bay_zips$ZCTA5CE10) %>% 
   mutate(
     DTE =
       as.Date(paste0(YEAR,"-",MONTH,"-1"))
     
   ) %>% 
  group_by(DTE, CUSTOMERCLASS) %>% 
  summarize(
    TOTALKBTU = 
      sum(
        TOTALKBTU, 
        na.rm = T
      ),
    TOTALCUSTOMERS =
      sum(
        TOTALCUSTOMERS,
        na.rm = T
      )
  ) 

bay_precovid <-
  pge_consump %>% 
  filter(ZIPCODE %in% bay_zips$ZCTA5CE10) %>% 
  filter(
     YEAR %in%
       c(2017, 2018, 2019), 
      MONTH %in%
        c(03, 04, 05, 06)
  ) %>%
  group_by(ZIPCODE, MONTH, CUSTOMERCLASS)%>% 
  filter( 
         CUSTOMERCLASS %in% 
           c("Elec- Residential","Elec- Commercial")
         )  %>% 
  summarize(
    TOTALKBTU = 
      sum(
        TOTALKBTU, 
        na.rm = T
      )/3,
      TOTALCUSTOMERS =
      sum(
        TOTALCUSTOMERS,
        na.rm = T
      )
  ) %>% 
  rename(
     PRECOVIDKBTU = TOTALKBTU
  )

bay_covid <-
  pge_consump %>% 
  filter(ZIPCODE %in% bay_zips$ZCTA5CE10)%>%
  filter(
     YEAR %in%
       c(2020),
     MONTH  %in%
       c(03, 04, 05, 06)
     
   ) %>% 
  group_by(ZIPCODE, MONTH, CUSTOMERCLASS)%>% 
  filter( 
         CUSTOMERCLASS %in% 
           c("Elec- Residential","Elec- Commercial")
         )  %>% 
  summarize(
    TOTALKBTU = 
      sum(
        TOTALKBTU, 
        na.rm = T
      ),
    TOTALCUSTOMERS =
      sum(
        TOTALCUSTOMERS,
        na.rm = T
      )
  ) %>% 
  rename(
     COVIDKBTU = TOTALKBTU
  ) 

final<- left_join(bay_precovid, bay_covid, by = c("ZIPCODE", "MONTH", "CUSTOMERCLASS") , copy = FALSE)%>% 
  mutate(
    NETKBTU= 
      COVIDKBTU - PRECOVIDKBTU
  )

Graph Observations

Graphing the total kBTU consumption showed a pattern of energy consumption that generally followed trends of previous years. There’s a general increase of residential consumption and decrease of commercial consumption, but the difference seems to mitigate overall consumption levels.

pge_chart <-
  bay_consump %>% 
  ggplot() +
  geom_bar(
    aes(
      x = DTE %>% factor(),
      y = TOTALKBTU,
      fill = CUSTOMERCLASS
    ),
    stat = "identity",
    position = "stack"
  ) +
  labs(
    x = "Date",
    y = "kBTU",
    title = "PG&E Bay Area Monthly Electricity & Gas Usage, 2017-20",
    fill = "Utility Type"
  )

pge_chart

pge_chart %>% ggplotly()

# Map Observations & Reasoning

For my maps examining difference in electricity consumption, I focused specifically on the month of April.

—I chose April because this is when the strictest shelter-in-place rules were enacted across the Bay Area. These rules did not go into effect until at least 10 days into March, and as months passed, each county became less restrictive. For these reasons, I thought examining April would show the most explicit differences due to COVID.

—I averaged electricity consumption (by consumer type) in the month of April from 2017-2019, and I subtracted that from the electricity consumption in April 2020. I chose to split residential and commercial consumption into two maps, as the trends seemed to go in opposite directions for the two types.

—Residential consumption generally increased , as more individuals were forced inside. However, these increases are more significant in zipcodes around the mouth of the Bay: San Francisco, Palo Alto, San Jose, and less so in zipcodes further out. I would attribute that to both the number of residents and the affluence of these zipcodes (more affluent locations generally contain less essential workers and can afford to stay home for longer.)

—Commercial consumption generally declined, as more businesses were forced to shutter.

—An important caveat is that this data does not examine electricity consumption on a more individual level, just the overall difference. For this reason, it is difficult to tell whether the differences in consumption are due to more people using electricity or a few more affluent people using a lot more electricity in certain zipcodes.

final_res_april <-
  final %>% 
  filter(
    CUSTOMERCLASS == "Elec- Residential",
    MONTH == 4)%>% 
  mutate(
    ZIPCODE = ZIPCODE %>% as.character()
  ) %>% 
  group_by(ZIPCODE) %>% 
  right_join(
    bay_zips %>% select(GEOID10),
    by = c("ZIPCODE" = "GEOID10")
  ) %>% 
  st_as_sf() %>% 
  st_transform(4326)

res_pal <- colorNumeric(
  palette = "Blues",
  domain = 
    final_res_april$NETKBTU
)

leaflet() %>% 
  addTiles() %>% 
  addPolygons(
    data = final_res_april,
    fillColor = ~res_pal(NETKBTU),
    color = "white",
    opacity = 0.5,
    fillOpacity = 0.5,
    weight = 1,
    label = ~paste0(
      round(NETKBTU), 
      " kBTU difference ",
      ZIPCODE
    ),
    highlightOptions = highlightOptions(
      weight = 2,
      opacity = 1
    )
  ) %>% 
  addLegend(
    data = final_res_april,
    pal = res_pal,
    values = ~NETKBTU,
    title = "Difference between<br> in Residential Electricity Usage  <br> b/w average of Apr 2017-2019 <br> &  Apr 2020, kBTUs <br>"
  )

final_com_april <-
  final %>% 
  filter(
    CUSTOMERCLASS == "Elec- Commercial",
    MONTH == 4)%>% 
  mutate(
    ZIPCODE = ZIPCODE %>% as.character()
  ) %>% 
  group_by(ZIPCODE) %>% 
  right_join(
    bay_zips %>% select(GEOID10),
    by = c("ZIPCODE" = "GEOID10")
  ) %>% 
  st_as_sf() %>% 
  st_transform(4326)

res_pal <- colorNumeric(
  palette = "Reds",
  domain = 
    final_com_april$NETKBTU
)

leaflet() %>% 
  addTiles() %>% 
  addPolygons(
    data = final_com_april,
    fillColor = ~res_pal(NETKBTU),
    color = "white",
    opacity = 0.5,
    fillOpacity = 0.5,
    weight = 1,
    label = ~paste0(
      round(NETKBTU), 
      " kBTU difference ",
      ZIPCODE
    ),
    highlightOptions = highlightOptions(
      weight = 2,
      opacity = 1
    )
  ) %>% 
  addLegend(
    data = final_com_april,
    pal = res_pal,
    values = ~NETKBTU,
    title = "Difference between<br> in Commercial Electricity Usage  <br> b/w average of Apr 2017-2019 <br> &  Apr 2020, kBTUs <br>"
  )